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Abstract 

In this paper we show how, under certain restrictions, the hydrodynamic equations for the freely 
evolving granular fluid fit within the framework of the time dependent Landau-Ginzburg (LG) 
models for critical and unstable fluids (e.g. spinodal decomposition). The granular fluid, which is 
usually modeled as a fluid of inelastic hard spheres (IHS), exhibits two instabilities: the spontaneous 
formation of vortices and of high density clusters. We suppress the clustering instability by imposing 
constraints on the system sizes, in order to illustrate how LG-equations can be derived for the order 
parameter, being the rate of deformation or shear rate tensor, which controls the formation of vortex 
patterns. From the shape of the energy functional we obtain the stationary patterns in the flow field. 
Quantitative predictions of this theory for the stationary states agree well with molecular dynamics 
simulations of a fluid of inelastic hard disks. 
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1 Introduction 



Granular matter Q consists of small or large macroscopic particles. When out of equilibrium, its 
dynamics is controlled by dissipative interactions, and distinguished in quasi-static flows or granular 
solids on the one hand, and rapid flows or granular fluids (2[ on the other hand. 

Typical realizations of granular solids are sand piles, avalanches, Saturn's rings, grain silos, stress 
distributions. Here particles remain essentially in contact, and the dynamics is controlled by gravity, 
friction and surface roughness. In this paper we concentrate on granular fluids. Typical examples 
are driven granular flows, such as Poisseuille flow vibrated beds |§, |5[ ^, or rapid flows with some 
form of continuous energy input |^]. Here the dynamics is controlled by inelastic binary collisions, 
separated by ballistic motion of the particles. The forces are of short range and repulsive, and the 
system is frequently modelled as a collection of smooth inelastic hard spheres (IHS) [|| of diameter a 
and mass m. Momentum is conserved during collisions, which makes the system a fluid, but energy 
is not conserved. In a collision, on average, a fraction e of the relative kinetic energy of the colliding 
pair is lost, where e is referred to as the degree of inelasticity. In the literature H, ||, e = 1 — a^, 
usually expressed in terms of the coefficient of restitution a. Its detailed definition does not concerns 
us here. 

Here we focus on the idealized limiting case of a freely evolving rapid flow without energy input 
and with nearly elastic collisions, and therefore slowly cooling. This system shows an interesting 
instability. When prepared in a spatially homogeneous equilibrium state, the system does not 
stay there, but slowly develops patterns, both in the flow field (vortices), and in the density field 
(clusters), the so called clustering instability. For the two-dimensional case the analogies with 
spinodal decomposition have already been pointed out in the literature pO| . 

The search for the proper macroscopic description of unstable granular fluids has been pursued 
by many authors [1-13,18,20,22-32]. Recently, two new points of view have been presented, namely 
by Ben-Naim et al. Q and by Soto et al. |l^. The flrst one needs to be discussed in more detail 
because macroscopic equations for granular fluids like the Burgers equation appear in that paper, 
as well as in the present one. These authors conjecture that the (vector) Burgers equation describes 
the flow velocity u(r, t) of a granular fluid, or at least of a dilute granular gas -the authors are not 



very explicit on this point 1 13 - on large space and time scales for arbitrary values of the inelasticity 
(0 < e < 1) in dimensions, i.e. 

= -u • Vu 4- I'V^u, (1) 

where v is the kinematic viscosity, and the solutions of interest have to satisfy V x u = 0. An 
interesting aside is that the rotation-free solution of the vector Burgers equation can be expressed 
as the gradient of a scalar fleld, u — V/i(r, t), where h{r, t) satisfies the equation, 

dth^vV^h+\\Vh\^. (2) 

This equation is the famous KPZ-equation, named after Khadar, Parisi and Zhang [Q. It describes 
the growth dynamics of solid surfaces, where h is the height function. 
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It is well known that the Burgers equation in the inviscid limit (i/ — > 0) describes the "sticky 
dust" or adhesion model of perfectly inelastic point particles with e — 1 [ p^ , [l7t . The conjectured 
validity of the Burgers equation for the inelastic hard sphere fluid implies a universality hypothesis, 
i.e. the large (r, i)— behavior of a granular gas is in the universality class of sticky dust, independent 
of its inelasticity e. Several implications of this conjecture in dimensions d > 2 have been criticized 
m Ref H, For the one-dimensional case the authors present rather convincing evidence based 
on molecular dynamics (MD) simulations using TV — 10^ particles, and on scaling arguments. In 
addition, Boldyrev ||ig|] has studied the one-dimensional randomly driven Burgers equation with a 
pressure term — (Vp)/p with p cx included. He has shown that properties, like the structure 
factor or pair correlation function, are not affected by the pressure term if a < 2. On account 
of this argument, the Vp-term in the one-dimensional Navier-Stokes equation can be neglected 
on the largest (r,t)~ scales for a dilute gas in the inviscid limit, making the Burgers equation 
the appropriate macroscopic equation. Note that in the freely cooling granular fluid the viscosity 
decreases as Vt as T — 0. In higher dimensions the systems used in MD simulations {N — 5 x 10^ 
in two dimensions and N — 10^ in three dimensions |Q) are too small to draw any conclusions 
on large scale behavior. Here analytic or scaling arguments to support the conjecture are lacking. 

Although equations of the form (|l|) will frequently appear in the present paper, we emphasize that 
Ref. [0 refers only to the largest possible scales (where the thermodynamic limit has been taken). 
In this paper we explicitly restrict ourselves to small systems in order to suppress the clustering 
instability. This has been done to simplify the problem. Moreover, Ref. 11 1 only refers to dilute 
granular gases, whereas this paper covers also liquid densities. 

A second new development, which is somewhat similar to ours, is given by Soto et al. [ p^ . 
These authors also study granular fluids, contained in systems that are sufficiently small, such 
that the clustering instability is suppressed. In these small systems the growth of vortices is very 
slow. Consequently, the remaining hydrodynamic modes are enslaved by the slowly growing unstable 
vorticity modes and the amplitude of these vorticity modes remains very small. Under this condition, 
they obtained the amplitude equations for the slowest vorticity modes, which can be derived from a 
potential function. 

The question of interest in the present paper is: can the models for granular fluids be fitted into 
the generic classification of Landau-Ginzburg-type models, as given by Hohenberg and Halperin 
to describe critical dynamics and hydrodynamic instabilities? The goal of this article is to illustrate 
how, under certain restrictions, the standard nonlinear hydrodynamic equations for the IHS fluid 

^ can be cast into a Landau-Ginzburg-type equation of motion for the order parameter, which 
can be derived from an energy functional and, more specifically, to point out which terms in the orig- 
inal hydrodynamic equations are responsible for the quartic terms in the Landau-Ginzburg energy 
functional. 

The plan of the paper is as follows. In Sec. 2, we start with the hydrodynamic equations. The 
decay of the total energy at short times and the results of a linear stability analysis are briefly 
reviewed. In Sec. 3, we introduce an assumption of incompressible flows under certain restrictions on 



3 



system size or time regime. Then, under these assumptions, the hydrodynamic equations are reduced 
to a closed equation for a scaled flow field. It is shown in Sec. 4 that this equation for a scaled flow 
field can be cast into the form of a time-dependent Landau-Ginzburg equation for an appropriate 
order parameter. The shape of the energy functional is discussed and possible stationary solutions 
are presented. Finally, in Sec. 5, we make a quantitative comparison of the theoretical predictions 
at large times with molecular dynamics simulations of inelastic hard disks. We end with some 
conclusions in Sec. 6. 



2 Dynamic Equations and Instabilities 

The starting point are the hydrodynamic equations. They are not only needed here to recapitulate 
our present theoretical understanding of this system, but also to formulate the new extensions to be 
discussed in this paper. 

The macroscopic time evolution of the IHS fluid on large spatial and temporal scales ||, ^ can 
be described by the nonlinear hydrodynamic equations for the local density n(r,t), the local flow 
field u(r,t) and the local temperature T(r,t), supplemented with a sink term F accounting for the 
energy loss through inelastic collisions, i.e. 

Dtn = — uV • u, 

Dtu = - — Vp + 2//V-D, 

mn 

DtT = -^V-u + 6TV^r + 26^D : D-F. (3) 
an 

In this paper the inelasticity is always assumed to be small. For later convenience the macroscopic 
equations are given for a d-dimensional systems. Here = 9* + u • V is the time derivative in 
a comoving frame. The local energy density of the IHS fluid is e = ^mnv? + ^nT, and p is the 
pressure. The shear rate I?q/3 is the symmetrized dyadic, {VqW^}, which is also made traceless. 
The coefficient 6t = 2k/ dn is proportional to the heat conductivity k, and h±^ — 2mv/d to the 
shear viscosity v. For simplicity of presentation the bulk viscosity has been set equal to zero, and 
the transport coefficients v and k, which depend on the local density and temperature, are taken at 
some fixed reference values, n{t) and T{t), to be specified later. The four terms in the energy balance 
equation account for work done by the pressure, for heat conduction, for nonlinear viscous heating 
and collisional dissipation. Gradients in the flow velocity considerably slow down the coUisional 
cooling process through viscous heating. 

On the basis of kinetic theory one can derive that the rate of collisional energy loss, F ~ 2^oU!T, 
is proportional to the collision frequency uj multiplied by the fraction of energy eT lost per collision 
where 70 — e/2d = (1 — a'^)/2d. In general, the collision frequency uj{T) is proportional 
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to the root mean square velocity vq = y^T/TTi, and its explicit form for hard sphere fluids can be 
found in Refs. |2^. When the system is prepared initially in a homogeneous equilibrium state, 
it evolves at short times in a spatially homogeneous cooling state (HCS) with a time dependent 
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temperature. Combination of Eqs. (^) and the expression for F given above, yields dtT 70^^/^. 

This leads to Haff 's homogeneous cooling law |^ for the mean energy per particle in the IHS fluid, 

m = lm^^^^^^=E,eM-2,or), (4) 

which is needed for later comparison. Here Eq — {d/2)To is the initial energy, and Iq — 1/ujq with 
luq = Lu{To) is the mean free time in the initial state. In general the number of collisions per particle 
T(t) in a time t is defined through dr = uj{T{t))dt. In the HCS integration of this relation yields 
exp(7or) = (1 +-foUJot). 

However, this state is unstable against spatial fluctuations in density n{r, t), temperature r(r, t) 
and flow velocity u(r, t). The present theoretical understanding of these instabilities 29, 2^, 

^ is essentially based on a linear stability analysis of the hydrodynamic fluctuations in the density, 
Sn = n ~ n, temperature, ST = T — T, and flow velocity u. This is done by using the rescaled 
Fourier modes 6ny^, 5T\^ = 5T\^/T and Jut ^ Uk/vT", where an overline denotes a spatial average, 
a = /y rfr a(r), and T is the global granular temperature. Fourier transforms are defined as 

/k = /y dr e^*'' ''/(r). The rescaled eigenmodes are described by 5a\^{T) — eyiY>{z\T)5a\^{Q). The 
exponential growth rates of unstable {zx{k) > 0) and stable {z\{k) < 0) modes are shown in Fig. 1 
as a function of the wave number k. 

This figure shows that the transverse flow field uj_k or shear mode (A =J-) with a wave number 
fc < fc^ is unstable, and develops vortices. On the other hand density fluctuations couple weakly, 
in order 0{k), to the heat modes (A = H), and znik) in Fig.l shows that these fluctuations are 
unstable in the range k < k*^^ and linearly stable in the range k > k*^ , i.e. remain at thermal noise 
level. The stability thresholds k*^ and k*j^ are defined as the the root of z\{k) — for A = 
and are marked as black dots in the figure. The figure also shows that the growth rate z^{k) for 
the vorticity mode is much larger than the growth rate for the heat mode znik), which couples to 
the density fluctuations. This explains why vortices appear long before the density clusters start to 
appear. 

An intuitive explanation for the appearance and growth of large scale vortices is that a binary 
collision destroys a fraction of the kinetic energy of relative motion of the colliding pair. The 
cumulative effect of many successive collisions is that they make the particles move locally in a more 
parallel and coherent fashion. This creates locally patches of vorticity. These patches grow in size 
by selective suppression (stronger damping) of short wavelength modes . 

An intuitive explanation of the appearance of high density clusters goes as follows |^ . Suppose 
one prepares the system initially in a spatially homogeneous equilibrium state, and there occurs 
locally a spontaneous negative pressure fluctuation ('depression'). The resulting particle flow from 
the surroundings tries to compensate for the local depression. This creates an excess local density, 
which increases the collision frequency, and in turn decreases the temperature. This creates again 
a depression, and the process keeps repeating itself, thus creating cold dense clusters, surrounded 
by a hot dilute gas. Moreover, these arguments also suggest that the pressure fluctuations and 
the corresponding gradients remain substantially smaller than those in density and temperature. 
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Figure 1: Dispersion relations Z\{k) (from right to left) for the shear mode (A =-L), heat mode (A = H) and 
the sound modes in the IHS fluid for an area fraction of (j) = 0.4 and a restitution coefficient a = 0.85. The 
stars mark the location of the minimum wave vector allowed in the system, kg = 2n/L, for the number of 
particles indicated, and the black dots mark the location of the threshold values k* ~ with corresponding 
correlation lengths where a — {\\,H,1.). 

In fact it would be very interesting if one could test this suggestion by measuring the pressure 
locally, i.e. in cells of the typical size of a density cluster, through molecular dynamics simulations, 
or Direct Monte Carlo Simulations (DSMC) of the Boltzmann or Enskog equation, or by the Lattice 
Boltzmann method. 

The shape of the dispersion relations for the growth rates also explains the suppression of insta- 
bilities through a reduction of the system size |2|, |2^. Let kQ{N) ~ 2tt/L be the smallest 
wave number, allowed in a system of linear extent L at fixed density and with periodic boundary 
conditions. In systems with k*^ < ko{N) < k]_, there are growing vorticity modes, but all density 
fluctuations are stable according to a linear stability analysis. In systems with ko{N) < k*fj the fluc- 
tuations in the density and in the flow field are unstable However, systems with k^){N) > k\ 
do not show any instability. The smallest allowed wave numbers ko{N) for different numbers of 
particles, N = 50, 100, 200, . . . , 1000 at fixed packing fraction (f>, are shown as stars in Fig.l. This 
information on small systems will be used in later sections. 

We finally remark that several authors have also studied nonlinear terms in the macroscopic 
equations for granular fluids, such as the viscous heating term ^ and the nonlinear 

convective term [|ll| . The inclusion of the combined effects of viscous heating and coUisional cooling 
is essentially the new mechanism driving the dynamics of dissipative granular fluids in the time 
regime, directly following the homogeneous cooling state. 

In an unpublished report one of us has developed a systematic expansion method for solving 
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the nonlinear hydrodynamic equations for granular fluids. The method is based on the separation 
of time scales of the relatively slow process of vorticity diffusion and the relatively fast process of 
heat conduction for all fc-dependent excitations, allowed in the system (fc > ko{N)). This picture is 
consistent with the relative sizes of the growth rates z±{k) — 70 and znik) — 70, as shown in Fig. 1. 

The separation of the time scales permits a systematic approximation scheme, based on the 
Chapman- Enskog method or multi-time scale expansion method Q . In zeroth order this method 
yields a closed equation of motion for the rescaled flow field u, as well as an equation for the global 
granular temperature T{t). In first approximation (see Ref. j32|) this method yields equations for the 
density and temperature fluctuations, 5n and 6T, and a higher order correction to the Navier-Stokes 
equation. 

In the present paper we will not explain the rather technical calculations based on the multi-time 
scale method, but we will try to elucidate the essential physics by deriving the zeroth order results 
in a more intuitive fashion. The systematic derivation will be published elsewhere. 



3 Incompressible flows 

As discussed in the previous section, instabilities and pattern formation occur in two different local 
fields, u and n, and on two different time scales, namely first patches of vorticity appear, and only 
much later density clusters appear. As explained above, the appearance of density clusters can even 
be further delayed, or all together suppressed by decreasing the system size ||2^, |8[ ^ These 
observations suggest to analyze first the simplest case where fluctuations in density and temperature 
remain small. This would happen in the time regime following the unstable homogeneous cooling 
state, and possibly in the full time regime for sufSciently small systems, as suggested by the linear 
stability analysis. Of course the stability of solutions on the largest time scales is determined by 
the full nonlinear equations. The conditions for small n— and T— fluctuations might be realized in 
incompressible flows, where V • u = 0. Then the density remains constant in the comoving frame, 
and the temperature balance equation greatly simplifies. As is well known from standard fluid 
dynamics and from the theory of turbulence | ]35| , ^ , flows of elastic fluids are quite incompressible. 
This implies, 

V • u = or = k • Uk = 0, (5) 

where W||k is the longitudinal Fourier mode. Moreover, MD simulations and the theory of hydrody- 
namic fluctuations in granular flows |2j, ^ show that the incompressibility assumption, U||k — 0, 
remains valid down to very small wave numbers, satisfying fc^n ^ 1 , and ultimately breaks down at 
the largest wavelengths, where ^|| ^ I/70 is the largest intrinsic dynamic correlation length in IHS 
fluid. It satisfles the inequality, ^|| ^ = {v/uj^i^Y/'^ for nearly elastic systems. Both correlation 
lengths are indicated in Fig.l and deflned more explicitly in j25|] . 

Therefore, as a zeroth approximation to our nonlinear theory, we make the incompressibility 



assumption, V • u = 0, following Refs. |2J, |25|, and the equations of motion become. 
An = 0, 
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Dtu = Vp + j/V^u, 

mn 

DtT = hrV^T + 6_L [Vu + (Vu)^] : Vu - (6) 

Consequently n(r, t) is constant in a comoving frame, and the set of nonlinear equations (^) can 
not describe the growth of inhomogeneities in the density field, but supposedly describe the time 
evolution of the system as long as the density fluctuations are small. 

As a consequence of the incompressibility assumption the local density stays constant in time, 
and neither depressions, nor hot and cold regions can develop. Therefore, the pressure gradient in 

remains negligibly small as well, and the Navier-Stokes equation becomes, 

dtM -u • Vu + j/V^u, (7) 

where the flow velocity u = is purely rotational with V • u = 0. We also note that the diver- 
gence of Eq.(0), in combination with V • u = 0, reduces to {V aUp){y pUa) = at all times. As 
already discussed in the introduction, the Burgers equation for a rotation-free flow field has been 
conjectured as the large scale macroscopic equation for a d-dimensional dilute granular gas of 
arbitrary inelasticity e. 

Next we consider the temperature balance equation, which involves two processes: the diffusion 
process of heat conduction, where Fourier modes Tk decay with a rate brk^ , and the global process of 
coUisional cooling and nonlinear viscous heating, describing the decay of Tk for A: — > 0, or equivalently 
of T[t) = (l/V) Jy drT{r,t), referred to as global temperature. 

To split off the behavior of f(t) from Eq.(||) we take the spatial average of the T~equation, 
yielding for the global temperature, 

dtT^bA_W^~2joLof, (8) 

where overlines denote spatial averages. We have introduced the notation |Ap = J^ap l^a/sP, for 
a second rank tensor A. The transport coefficient b±^ and collision frequency to are functions of the 
spatially average density n — N/V and temperature T{t) (see 'reference values' below Eq. (||)). This 
is allowed as long as the local fluctuations Sn = n — n and ST — T ~ T are not getting too large 
through nonlinear effects. 

The new evolution equations and (|^) contain the time dependent coefficients b±, uj and 
which are proportional to uo(0 = {2T{t) /niY/'^ . Therefore, it is convenient to introduce the scaled 
field u — u/uo, and the scaled time r, defined as dr = u){T)dt. The final macroscopic evolution 
equations then become, 

9.1nT = ^I?i|Vii|2-27o, 

i9^u + /oU-Vu = Pj_V2u- i(a^lnT)u 

2 

= 7ou-f X>_lV2£i- -X>_l|Vu|2u, (9) 

d 

where = vq/oj is the (time-independent) mean free path. The last equation is a closed equation for 
u, and the physically consistent solutions need to obey the relation V-u = 0. The global temperature 
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is slaved by the flow field. The rescaled vorticity diffusion coefficient, defined as 'D± — vjuj, is 
independent of time. The first term on the right hand side of the equation for u accounts for the 
instability, the second for the vorticity diffusion and the last one for the saturation effects, caused by 
the nonlinear viscous heating. It slows down the growth of unstable k-modes, and may ultimately 
lead to a steady state for u. For large times the above equations and (|8|) arc no longer valid, and 
nonlinear effects will induce density inhomogeneities, even in small systems with k*^ < kQ{N) < k^. 



4 Spontaneous symmetry breaking 

In this section we will drop the nonlinear convective term u • Vu, but at the end of our analysis we 
admit out of all possible solutions only those that satisfy Eq. with the convective term included. 

The final equation for the rescaled u— field has the form of the Landau- Ginzburg equation of 
motion for a non-conserved order parameter. This can be made more explicit by introducing the 
order parameter S = Vu, and applying V to the u-equation in (^) with the result, 

drS - (7o+2?±V2)S-^I?i|SpS 

= -vsn[S]/6s\ (10) 

where the last line contains the functional derivative of the energy functional 7i[S], defined as 

n[S] = -Ij.jsi^+lv^WW+^'Di- {WY- (11) 

In our further considerations it is more convenient to use Fourier modes. Moreover Uk = u^k, 
because of the assumption of incompressibility, and Sk = kuj^k- Then we obtain the equation of 
motion, 

9,Sk = -V^Sn[S]/5Sl 



7o-I?^fc^-^ENp|sk (12) 



with an energy functional. 




^[SH^E(-^o+^^^')iSkp+^i— >:iSkn , (13) 

where the wave number k = does not contribute. 

The energy in (^l]) and ( p^ ) resembles a Landau free energy form for a tensorial order parameter 
S = Vu with trS = V • u = 0, which is in fact the shear rate or rate of deformation tensor. It has 
a quartic term S"^, and a quadratic term S'^ with a coefficient that vanishes at a critical threshold 
k'^ = (7o/2?±)^^^. It differs from the standard Landau free energy in that the quartic term contains 
summations over two totally independent wave numbers. 

These results are very interesting. If the energy functional has a minimum, then there is a fixed 
point solution, Sk(oo), that is approached for large times. They are found by setting the right hand 
side of ( p^ ) equal to zero, i.e., 
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Figure 2: Landscape plot of the energy Ti.[S] in the Sko ^subspace. 



(14) 



q 



We will show that depending on the parameter, ko — 2tt/L, being above or below the threshold 
value k*^_ = \/jq/ D± , the fixed point value of the order parameter, {Sk(oo)}, is vanishing or non- 
vanishing. A stable fixed point with some non-vanishing Fourier components indicates that the 
system approaches a non-equilibrium steady state with a stationary pattern in the fiow field, and 
spontaneous symmetry breaking has occurred. 

Consider the right hand side of (|l^), and observe that the expression between curly brackets is 
necessarily negative for k > k*^, and the Fourier mode Sk decays to zero. If the smallest possible wave 
number satisfies, fco > A:^, all Sk decay to zero, and there is no spontaneous symmetry breaking. 
The system remains spatially homogeneous. However, if fco < fc^, then there exists the possibility 
that the expression inside brackets in Eq.([l^) vanishes for a non- vanishing value of Sk(oo), i.e. there 
is an extremum determined by the condition. 



If the extremum is a saddle point then there are directions of exponentially growing solutions. 
However, if fc < fc^, then the right hand side of (|l5|) may become zero and even positive. There 
is the possibility of stationary and of exponentially growing solutions. The fixed point {Sk(oo) — 
for any k } is a saddle point. One can also show |32| that all fixed point solutions with non- 
vanishing Sk for any |k| ^ fco are saddle points, and that the only non-vanishing solution {Sk 7^ 
if |k| = fco, Sk = otherwise} is a stable fixed point with an infinitely degenerate minimum, given 
symbolically by the Mexican hat shape as illustrated in Fig. 3. The condition (|l^) for uj^k with 




(15) 



q 
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|k| = feo = 27r/L in d-dimensions is then 

1 d 

E l"koj' - 4(70 - -D^kD/V^kl (16) 

a— 1 

where kpa — k^ea and Gq, is a unit vector in the direction a {a = {1, 2, • • • d}). As the solutions of 
these equations have to satisfy Eq. (|^), the Fourier components can be written as 

fiko„ = e;3A„^e'«"^ (17) 

where the phases Oap and amphtudes A^f) are 2d{d — 1) real numbers. On account of ( p^ the 
amplitudes satisfy the relations 

^0 ^ E E ^'/3 = <To - T^^kD/Vj^kl (18) 

The stable stationary solutions in real space are then 

iia{v) =^^Aai3epcos{kara+ Oap) ■ (19) 

Out of this set of solutions we select those that satisfy the full nonlinear Eq. (||) with the convective 
term included, i.e. we determine the d{d — 1) amplitudes Aap such that the relation, uq • Vuq — 0, 
is satisfied. This yields d conditions. For the two-dimensional case only two fixed point solutions 
remain, instead of infinitely many degenerate minima, i.e. 

Uo(a;) = eyAQCOs{kf)X + 9x) 

UQ{y) = exAocoa{koy + 9y). (20) 

The symmetry of the steady state is spontaneously broken, and spontaneous fluctuations in the 
initial state determine which of these two minima will be reached. 

In summary, the equation, describing the growth dynamics of vortices in granular fluids, without 
the convective term is described by a time dependent Landau-Ginzburg equation for a non-conserved 
order parameter, S — Vu, derived from an energy functional with a continuous set of degenerate 
minima, having the shape of a Mexican hat. The last observation may suggest that unstable granular 
fluids have some resemblance to Model H which has a similar energy surface in the neighborhood 
of its stable flx points. However, this is not the case. Addition of the nonlinear convective term 
to Eq.(p^) selects out of this infinite set of minima only subsets of admissible solutions. In two 
dimensions only two distinct minima survive. Therefore, the unstable two-dimensional IHS fiuid 
has a greater resemblance to spinodal decomposition for a non-conserved order parameter, which 
belongs to the universality class of Model A [|T] . 

It should be noted that the complete solution of Eq. (^) with the convective term included may 
have a larger set of physically acceptable solutions. However, we have not been able to find any. 
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5 MD- Simulations 



In this section we present results obtained from the zeroth approximation and compare the theoretical 
predictions with the results of computer simulations of small systems in two dimensions. At the end 
of this section, some results of the first order approximation in Ref. are presented without 
derivation and compared with computer simulations. 

A snapshot of a typical configuration for a small system with < ko < k\ at large times 
T 3> Tcr is shown in Fig. ^. A shear flow with variation of the u^;— component in the y— direction 
is observed. Individual w^:— components of the particle velocities are plotted in Fig. ^ versus their 

coordinate. A fit to a sinusoidal curve (solid line) shows that the solution described in (^0|) is 
realized. 




Figure 3: Left: Configuration with A'^ = 400 at r = 600. A shear flow is observed in agreement with the 
solution of the nonlinear equations given in Eq. ([2o[). Right: Numerical data for this and subsequent plots. 
Last three entries are at r = 600. Units are chosen such that initial temperature To = 1, mass m = 1, and 
sphere diameter cr = 1. 



We assume that the number of collisions t in the simulations of Figs. 3 and 4 is sufficiently large, 



so that the components Uk are very close to their fixed point values. Then, the relation, Pj^jVup = 
\d{'^Q — I?j_fco), follows from ( p^ ) and (^, and c^t- Inf ~ —21)^^1. The global temperature for large 
times, i.e. t 3> = becomes, 

T'~Teexp(-2X'_Lfco^)- (21) 

Equation ( pT[ ) gives the temperature as a function of t. The integration constant cannot be de- 
termined from the theory, but a fit of ( pT]) to the simulation data in Fig. 5 for r > 200 yields 
Tg ~ 3.87 X 10"*. We also note that the mode coupling theory of Ref. |30 gives the same decay rate 
for the total energy, when this theory is applied to a small system. In that case the Fourier sum 
{1/V) can not be replaced by J (iq/(27r)''. On the contrary, the sum is essentially given by its 
first term only. 
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Figure 4: Profile of shear flow in Fig-^. Dots represent tlie instantaneous velocities of the particles; the solid 
line is a flt to the sinusoidal function, described in Eq. (po[), with amplitude Ag^^/vo ~ 2.78 and theoretical 
prediction Ao — 2.51. 



The relation between r and t is defined through dr = lu(T) dt, where the collision frequency to is 
proportional to the square root of T, i.e. 



/rr\ / ^ 

— = Lu{T) = iUo\ — 

dt V ^0 



T^f exp [-V^k^r] , 
Jo 



(22) 



where tuo = ui{To) is the collision frequency at the initial time. As oj can be measured directly 
as a function of t in MD simulations, it yields an independent determination of with the result 
Te ~ 3.68 X 10""*. Integration of (||) yields, 



(23) 



valid for {t — te) large. After eliminating r from Eq. ( pT| ) and (^3|) we obtain the global temperature 
for asymptotically large time t as, 

To 1 



(24) 



(L^02?±fc2)2 (t-te)^' 

where visual inspection of Fig. 6 shows that te — 10^. The temperature shows algebraic decay with 
the same exponent as in Haff's law (||), but the prefactor in Half's law does neither depend on the 
system size, nor on the dimensionality. In ( p^ it is proportional to L"^ ^ N^/'^ , whereas the prefactor 
in Haff's law is independent of the system size (see Fig. 6). 

Once we have T, we can calculate the energy per particle, defined in [BOl as. 



E 



-n 



-nT 





-TTT d 


= T 






" + 2 



(25) 
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Figure 5: Energy per particle and temperature as a function of number of collisions per particle r. At 
r — 600 one finds {E/T)^^^^ ~ 5.38 with a theoretical value 7o/I'±fco — 4.14 



where local density variations are supposedly small and the relation u = vqU has been used. Com- 
bination of this result with ( p^ ) yields, 

E~-^f. (26) 

This shows that at long times, the energy and temperature are proportional. Therefore, energy 
decays exponentially when expressed in terms of r (see (^)), with the same decay rate as the global 
temperature. This property is observed by MD simulations, as shown in Fig. |^. From ( |26| ) and ( p^ ) 
the decay of the energy per particle in terms of t in the long time limit is given by 

Hence, if we vary the number of particles N , at fixed packing fraction the amplitude of t^^- 
decay of the energy is proportional to iV'^ in two-dimensional systems, or to N^^"^ in d-dimensional 
systems. Comparison of Eq. (^7|) with simulations is shown in Fig. |^. It would also be interesting 
to compare the coefficient B in ( |27| ) with the simulation results of Ref . |lj] for the IHS fluid which 
have been performed in five and six dimensions. 

Finally, a systematic expansion to be presented in p2[ , shows that the stationary solution of 
local temperature and local density are given by 

6f^^-l -^cos[2(fcoy + e,)] 
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Figure 6: Left: Simulation results of the behavior of E{t) as a function of t for TV particles. At large times 
{t > 10" for <j} = 0.4 and a = 0.85), a power law tail ~ is observed. The coefficient B, defined in 
(p?!), is on average 2.3 x 10~^ with a theoretical value B ~ 1.2 x 10~^. Right: the amplitude of the tail is 
proportional to A^"^. Here the solid line is the result of the theory, Eq. ([27|), while the dots are simulation 
results extracted from the left panel. 

n(y) ^T^o'^fe) 
Sh=^-l ~ co8[2{koy + 0y)], (28) 



2duj 1 



where n — N/V, and 9y is the same phase factor as given in (^0|). It means that the density and 
temperature inhomogeneities show a period which is half the period of the shear flow profile. The 
relation between the spatial periods has already been observed in Ref. |^ in MD simulations of a fluid 
of two-dimensional hard disks, and in Ref. |^ in Direct Monte Carlo simulations of the Boltzmann 
equation for an IHS gas. Figure |^ shows the observed density and temperature inhomogeneities. A 
fit to a sinusoidal curve supports the temperature and density profiles given by (|2^). 

As shown in Eq.(|2^), the temperature and density inhomogeneities are in opposite phase, imply- 
ing that dense regions are cold, and the dilute regions hot. The amplitudes are such that the overall 
pressure is constant, as we have assumed in the course of the paper. 



6 Conclusions 

Under the restrictions imposed in our derivations, as discussed in Sec. 3, we have shown that the un- 
stable dynamics and formation of vortex patterns in the fiow field of a freely evolving fluid of inelastic 
hard spheres (IHS) can be cast in the form of a time-dependent Landau-Ginzburg-type model for a 
non-conserved order parameter. In the two-dimensional case (but not for d > 3) the growth of the 
vortex pattern in the IHS fluid is qualitatively similar to spinodal decomposition for a non-conserved 
scalar order parameter, referred to as model A in the Hohenberg-Halperin classiflcation pi[. The 
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Figure 7: Density (squares) and temperature (circles) inhomogeneities, Sh and ST at r = 600. The solid 
lines correspond to nonlinear fits to sinusoidal functions of half the period in ([28[). The simulated ratio of 
5T- to 5n amplitudes is here 2.0, and the theoretical prediction is 2.3. 

analogy between unstable IHS fluid in two-dimensions and spinodal decomposition has already been 
pointed out in Ref . |24) . A difference between our model and model A is that the energy functional 
( prj ) contains a quartic term with summations over two independent wave numbers. This implies 
a non-local interaction of the order parameter S = Vu. The non-local interaction is caused by the 
scaled field u = u/uq ^ u/Vt. Because the global temperature T is determined by S in all space 
(see Eq.(|[)), the evolution of a local order parameter S is affected by S at any other points in space 
through T. 

We have shown that nonlinear viscous heating gives rise to a quartic term in the energy functional, 
which is responsible for saturation of unstable vorticity modes. Unstable vorticity modes initially 
grow as predicted by the linear stability analysis, but eventually saturate because of nonlinear viscous 
heating. The influence of nonlinear viscous heating on the formation of temperature inhomogeneities 
and clustering has been pointed out by Goldhirsch and Zanetti |^, and investigated in more detail 
by Brey et al. by comparing the results of a hydrodynamic model with viscous heating, using 
direct Monte Carlo simulation of the Boltzmann equation. 

The theoretical predictions of our zeroth order theory for the flow field and the decay of the 
energy are in quantitative agreement with MD simulations of small systems as shown in Sec. 5. They 
support the intuitive arguments used in deriving the zeroth order description, presented in Sec. 3. 
Moreover, the results presented here can be obtained as the lowest order approximation in a multi- 
time scale expansion |]32| . This theory makes no assumption of incompressibility. There are some 
properties of those small systems that can not be described by the zeroth order theory. Density and 
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temperature inhomogeneities are only found in the next order of approximation. 

It is worth mentioning that results of our zeroth and first order approximations arc consistent 
with the results of a nonlinear analysis by Soto et al. ||T^ of systems which are close to the stability 
thresholds k]_ . If the smallest wave number of the system fco is slightly smaller than k]_ , vorticity 
modes with the wave number fcp grow so slowly that the rest of hydrodynamic modes are enslaved by 
these vorticity modes. Besides, the amplitudes of these vorticity modes at long times are expected 
to saturate because of nonlinear effects and remain small. On the basis of this consideration, they 
obtained amplitude equations for the vorticity modes with /cq. The inhomogeneous density and 
temperature are given as functions of the vorticity modes with ko. 

Finally, we emphasize that the periodic boundary conditions, with which the simulations were 
carried out, obviously play a crucial role in the formation of a shear flow profile and of the density 
and temperature profiles, presented in (|20| ) and ( p8|) respectively. Indeed, once the typical size of 
the vortex patterns becomes comparable to the length L of the system for t > Tcr — L'^/'D±, the 
artificial periodic boundary conditions start to affect the evolution of the system. Consequently, 
the long time regime r » Tcr described by the stationary solution is of less physical interest than 
the regime of unstable growth r ^ Tcr , where the typical size of the vortex patterns remains small 
compared to the length L of the system. A thermodynamically large system is always in the unstable 
growth regime. Unfortunately, the only analytic 'large' time results for the unstable growth regime 
of granular fluids in thermodynamically large systems, have been obtained from fluctuating (linear) 
hydrodynamic equations or from mode coupling theory |3^, but not from truly nonlinear theories. 
An exception is a dilute gas of inelastic point particles in one dimension , where strong evidence 
supports the conjecture that the large space-time behavior is described by the adhesion model and 
the Burgers equation. 
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